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ABSTRACT 


A  method  of  calculating  molecular  flow  properties  is  described, 
the  molecular  kinetics  (M.K.)  method,  which  is  distinctly  different 
from  the  Monte  Carlo  method.  It  has  the  advantage  of  permitting  the 
calculation  of  final  asymptotic  values  directly  from  initial  conditions 
once  the  transfer  probability  matrix  for  the  specific  geometry  is 
known.  Flow  histories  may  also  be  developed  by  a  variation  of  the 
method . 

The  general  theory  is  developed,  the  mathematical  operations 
are  described,  and  applications  are  shown  for  selected  two-  and  three- 
dimensional  configurations. 
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1.0  INTRODUCTION 

In  recent  years,  with  the  advent  of  requirements  for  space  environ¬ 
ment  simulation,  increased  attention  has  been  directed  to  the  behavior  of 
low  density  gas  flow  both  in  laboratory  experiments  and  theoretical  com¬ 
putations.  Especially  in  the  molecular  flow  regime,  computer  programs 
are  useful  tools  for  examining  interactions  between  a  gas  and  its  physical 
boundary  geometry,  the  walls  of  containers  and  ducts.  Two  distinct 
methods  are  in  use. 

The  first  method  is  the  popular  Monte  Carlo  technique  in  which  many 
individual  particles  are  followed  through  essentially  random  processes. 
Each  is  traced  over  a  complete  history  until  so  many  have  been  traced 
that  statistical  inferences  may  be  drawn,  extrapolating  to  the  behavior  of 
effectively  continuous  flow  of  large  numbers  of  particles. 

An  alternate  method,  of  prime  concern  here,  develops  from  radiant 
energy  transfer  analysis.  For  gases  in  the  free  molecule  regime, 
analogies  may  be  established  between  radiant  energy  flux  and  particle 
motion.  In  this  paper  we  call  this  method  "molecular  kinetics"  (M.  K. ). 
The  M.  K.  method  has  an  advantage  of  computer  speed  over  the  Monte 
Carlo  technique  as  long  as  the  geometrical  boundary  configuration  is  un¬ 
changed.  The  effort  required  to  change  configurations  is  about  equal  for 
the  two  methods. 

We  shall  refer  to  the  Monte  Carlo  method  only  for  comparison.  Both 
methods  share  a  number  of  assumptions  and  have  common  purposes,  that 
is,  to  determine  certain  gross  features  of  gas  boundary  interactions  at 
low  densities. 

The  following  fundamental  assumptions  are  used: 

1.  All  physical  surfaces  are  perfectly  diffusing;  Lambert's  Law 
holds . 

2.  No  intermolecular  collisions  occur. 

3.  An  ergodic  hypothesis  for  equilibrium  holds:  the  time  average 
velocity  of  one  molecule  is  equivalent  to  the  instantaneous  average 
velocity  of  a  large  number  of  molecules  so  that  a  single  velocity 
suffices  to  characterize  the  system. 

Similar  assumptions  occur  in  much  radiant  heat-transfer  theory, 
incorporated  in  assumptions  such  as  diffuse  reflection  and  emission, 
isotropic  non-dispersive  media,  and  the  use  of  geometrical  optics. 


Manuscript  received  May  1963. 
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2.0  MONTE  CARLO  METHOD 


A  geometrical  boundary  configuration  is  described,  initial  starting 
locations  (sources)  are  established  by  some  selection  rule,  and  the  line 
of  motion  of  a  particle  is  found  by  pseudo-random  methods  based  on  a 
statistical  distribution  model  such  as  Lambert's  Law.  The  equation  of 
the  line  of  motion  is  solved  simultaneously  with  equations  of  target  sur¬ 
faces  to  find  the  first  intersection.  This  point  is  treated  as  a  new  origin, 
and  a  new  path  line  is  generated.  The  tracing  stops  when  some  final  con¬ 
ditions  are  met,  such  as  "loss"  of  the  particle  back  to  the  source  or 
"capture"  by  a  sink.  Then  a  new  particle  is  generated  at  the  source. 

After  large  numbers  of  particles  have  been  traced,  one  may,  with 
caution,  begin  to  draw  statistical  inferences  about  the  behavior  of  effec¬ 
tively  continuous  streams  of  particles.  This  method  may  be  thought  of 
as  a  combination  of  discrete  particles  with  piecewise  continuous  geometry 
of  the  boundary.  The  molecular  kinetics  method,  in  contrast,  may  be 
thought  of  as  a  combination  of  continuous  flow  with  discrete  segmented 
geometry.  Both  methods  are  practical  approximations  to  the  real  prob¬ 
lem  in  which  both  boundary  and  flow  are  in  some  sense  continuous;  solu¬ 
tions  are  usually  not  available  for  this  problem. 


2.1  ERRORS 

The  precision  of  the  Monte  Carlo  technique  is  inherently  limited  only 
by  time  if  care  is  taken  to  reduce  the  possibility  of  systematic  bias.  The 
M.  K.  method  is  limited  in  precision  by  the  fundamental  assumption  of 
uniform  density  over  each  increment  of  area  at  each  step.  The  error  is  at 
its  worst  for  concave  corners  but  can  be  reduced  by  adopting  non-uniform 
area  increments  to  improve  the  approximations  where  the  densities  are 
changing  most  rapidly  with  position. 


2.2  COMPUTING  PROBLEMS 

The  Monte  Carlo  method  requires  time  principally  from  the  need  for 
large  numbers  of  runs  so  that  statistical  averages  have  validity.  The  ran¬ 
domness  of  the  process  is  never  overcome  completely.  The  computer 
program  acquires  bulk  when  the  boundary  of  the  system  is  to  any  degree 
complex,  and  it  requires  sophisticated  programming  to  avoid  long  search¬ 
ing  for  intercepts  of  lines  of  motion  with  the  boundary. 
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3.0  MOLECULAR  KINETICS  METHOD 


The  geometry  is  subdivided  into  discrete  elements  for  all  pairs  of 
which  "angle  factors"  (influence  or  form  factors)  may  be  calculated.  * 
These  factors  have  the  property  that  for  any  source  elements  facing  a 
closed  boundary  of  surface  elements,  all  factors  sum  to  unity.  They 
are  also  by  definition  a  distribution  function;  they  give  the  fraction  of 
radiation  leaving  the  source  element  and  intercepted  by  each  of  the  tar¬ 
get  elements.  They  may  thus  be  interpreted  as  transfer  probabilities, 
giving  the  probability  that  a  particle  from  a  source  element  will  strike 
a  target  element.  Finally,  they  are  functions  of  relative  orientation  be¬ 
tween  two  elements,  hence  functions  of  location  of  both.  We  shall  use 
the  name  "transfer  probability"  for  "angle  factor". 


The  geometrical  segments,  or  zones,  may  be  labeled  in  order; 
thus  Zi  is  the  i*h  zone  with  i  ranging  from  1  to  N,  a  finite  but  often 
large  number.  Then  Py  will  represent  the  probability  of  transfer  from 
zone  i  to  zone  j.  By  definition  of  Pjj  we  have 


,5.  P‘i  “  1 


(1) 


i.  e. ,  all  particles  leaving  arrive  somewhere.  The  Pjj  may  be  con¬ 
veniently  thought  of  as  a  square  matrix  of  size  n  x  n.  We  shall  use  the 
probability  matrix  to  distribute  particles  from  source  elements  i  to 
target  elements  j.  We  start  with  a  given  initial  distribution  of  particles 
in  one  or  more  zones.  Since  there  are  N  zones,  we  may  think  of  the 
original  particle  distributions  in  source  zones  as  the  N  components  of  a 
vector,  V^O),  i  =  1,  2,  .  .  . ,  N.  At  least  one  j8  non-zero.  By  this 

approach  we  easily  describe  any  original  particle  distribution  in  any  pat¬ 
tern  over  the  geometric  elements. 


We  have  implicitly  introduced  one  source  of  error  by  the  tacit  assump¬ 
tion  that  the  distribution  of  particles  within  a  zone  is  either  uniform  or 
has  some  assignable  mean  location.  It  is  plain  that  any  distribution  can  be 
approximated  as  closely  as  we  wish  by  making  N  sufficiently  large  so  that 
the  area  elements  approach  infinitesimal  size. 

The  original  distribution  Vj  is  transferred  as  a  pulse  to  all  target 
zones  by  the  simple  process  to  produce  the  next  distribution, 

vj(l>  =  J,  Pij  Vi(,)  (2) 


♦Jakob.  Heat  Transfer.  Vol.  2,  Chapter  31,  John  Wiley, 
New  York. 
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We  must  now  invoke  the  assumption  that  the  particles  arriving  at 
the  target  zone  are  uniformly  distributed  over  that  zone  if  we  are  to 
proceed. 


We  introduce  now  a  reflection  coefficient,  having  a  value  between 
zero  and  unity  and  dependent  on  that  zone,  which  we  label  A^.  If 
Ai  =  1,  the  zone  absorbs  no  particles,  but  if  Ai  =  0,  it  absorbs  all  the 
incident  particles.  By  this  means,  we  may  have  porous  walls  or  may 
introduce  holes  in  the  boundary.  Any  zone,  i,  for  which  Aj  =0,  is  a 
hole  from  which  incident  particles  do  not  return. 


We  may  now  calculate  one  step  further.  The  distribution  remain¬ 
ing  after  the  first  transfer  (Eq.  (2))  may  be  modified  by  the  reflection 
coefficient  so  that  the  resulting  distribution 


Aj  Vj<l)  =  A 


l  ,1.  Pii  Vi' 


(o) 


(3) 


remains  to  be  transferred  by  the  process 

Vj(,)  =  J,  Pij  Ai  v/1’ 


(4) 


The  repeated  application  of  the  above  steps  is  called  the  "bounce" 
approach  and  may  be  continued  indefinitely.  If  any  Aj  are  less  than 

the  original  particles  will  eventually  all  leak  out,  all  components 
becoming  vanishingly  small,  and  the  process  is  usually  terminated 
when  the  sum  (or  largest)  of  the  becomes  insignificant  relative  to 

the  original  total  distribution 


unity. 

Vi(m 


N 

2 

i  =  1 


Vi 


to) 


Several  quantities  of  interest  may  be  developed.  Suppose  that 
zone  K  is  a  hole-  that  is,  Afc  =  0.  Then  the  leakage  through  zone  K  at 
bounce  M  is  Vk'™'.  This  quantity  may  be  summed  for  all  bounces  to 
give  the  total  transmission  through  the  hole,  which  may  be  one  of  many 
or  part  of  a  larger  hole,  so  that  we  may  obtain  transmission  through  the 
hole  as  a  distribution  over  its  parts. 

Quantities  so  developed  approach  asymptotic  values  which  may  be  ob¬ 
tained  directly.  It  is  accepted  that  an  infinite  sum  over  a  history  of  one 
original  distribution  is  equivalent  to  a  sum  of  current  contributions  of  an 
infinity  of  like  distributions  originating  uniformly  in  the  past.  This 
assumption  allows  direct  development  of  asymptotic  values  and  implies 
a  steady  state  of  flow. 
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3.1  CALCULATIONS 

3.1.1  Asymptotic  Calculations 

We  modify  notation  slightly  so  that  Vq  represents  our  original  dis¬ 
tribution  vector,  P  is  the  matrix  of  elements  (Pij),  and  A  is  a  diagonal 
matrix  of  reflection  coefficients  (Ajj).  We  define  the  "hit"  vector  H  for 
molecules,  incident  on  the  zones  as  follows:  subscripts  denote  "bounce" 
number.  P'  is  the  transpose  of  P. 

H,  ■=  P'  V„  A  H,  =  V, 

H,  -  P'  A  H,  A  H,  =  V, 

H,  =  P'  A  H, 

Summing  over  all  bounces  for  Ht  =  H,  +  Ha  +  H,  +  •  •  •  • 

I  HT  =  P'  V,  +  P'  A  Ht 

where  I  is  the  identity  matrix  and  by  rearrangement, 

(I  -  P'  A)  Ht  =  P'  V„ 

Hence 

Hr  =  (I  -  P'  A)"‘  (P'  V.)  (5) 

for  the  asymptotic  value  of  total  hits.  Ht  is  a  vector,  and  its  components 
are  total  hits  per  zone. 

If  we  define  F  as  the  flux  of  particles  from  a  zone,  we  have  in  the 
same  way 

F.  =  V. 

F,  -  A  P'  V0 
F,  «  A  P'Fl( - 

or  I  F-r  -  A  P'  Ft  +  V,  +  A  P'  V.  (6) 

thus 

FT  •  (I  -  A  P')“'  (V,  ♦  A  P'  V.) 

The  hit  vector  H  is  related  to  pressure  distributions,  flux  vector  F  is 
related  to  the  density  in  the  space  near  the  boundary  within  the  system, 
and  their  difference  (H-F)  is  related  to  the  leakage  out  of  the  system  (i.  e.  , 
pumping  loads).  These  quantities  are  vectors  whose  elements  associate 
with  their  corresponding  zones. 
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3 .1.2  Computing  ProLiant 

In  the  molecular  kinetic  approach,  all  geometry  interactions  are  in¬ 
corporated  into  the  transfer  probabilities.  Herein  lies  the  bulk  of  com¬ 
puting  time.  These  probabilities  are  not  always  simple  to  compute,  as 
will  be  brought  out  in  later  applications.  But  once  computed  for  a  given 
configuration,  the  probability  matrix  need  not  be  repeated.  There  re- 
’mains  the  problem  of  computer  memory  size  since  the  P  matrix  is  to  be 
stored,  and  for  asymptotic  cases,  room  must  be  allowed  for  the  matrix 
inversion  and  storage.  Thus,  the  fineness  of  geometry  subdivision  is 
controlled  primarily  by  computer  storage.  This  problem  is  partly  solved 
by  block  programming,  logical  blocks  being  (1)  probability  calculations, 
(2)  preparation  of  matrix  for  inversion,  (3)  matrix  inversion,  and  (4)  final 
calculations  (probability  calculations  or  "bounce"  transfer  calculations 
repeated) . 

Block  programs  along  these  lines  have  been  prepared  for  limited 
general  cases  with  only  step  (1)  needed  for  specific  applications.  In 
addition,  other  programs  have  been  prepared  to  take  advantage  of  sym¬ 
metry  problems  and  matrix  partitioning  when  these  can  be  utilized  in 
special  applications. 

For  some  configurations,  we  have  in  the  molecular  kinetics  method 
a  much  faster  program  than  the  Monte  Carlo  approach  will  allow.  In  the 
following  sections  we  will  indicate  how  an  M.  K.  program  can  provide 
information  to  improve  "convergence"  of  a  Monte  Carlo  program. 


3.2  APPLICATIONS 

3.2.1  Two-Dinonslonol  Applications 

3 .2.1.1  Two-Oinansionai  Modal* 

Many  configurations  of  interest  may  be  described  approximately  by 
general  cylindrical  surfaces  generated  by  a  family  of  parallel  straight 
lines  of  infinite  length.  An  infinitely  long  strip  of  infinitesimal  width 
characterizes  a  surface  element.  If  all  surface  elements  are  of  this  kind, 
then  the  probabilities  for  particle  transfer  are  independent  of  location  along 
the  length  of  the  source  strip,  and  an  integration  may  be  performed  along 
the  target  strip  to  obtain  probability  of  transfer  from  source  strip  to  target 
strip.  The  configuration  may  be  represented  by  any  cross  section  perpen¬ 
dicular  to  the  generating  lines.  The  cross-section  plane  representation  of 
cylindrical  surfaces  is  what  is  usually  meant  by  a  two-dimensional  model. 
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Typical  of  current  uses  of  Monte  Carlo  and  M.  K.  programs  are 
calculations  related  to  the  baffle  configurations  used  with  cryopump 
panels.  These  low  temperature  pumping  surfaces  have  low  heat  con¬ 
duction,  and  it  is  desirable  that  molecules  arriving  at  the  cryosurface 
should  be  of  low  energy.  This  is  accomplished  by  insuring  that  a  large 
part  of  the  molecules  have  lost  energy  by  collisions  with  low  tempera¬ 
ture  (but  not  pumping)  baffle  surfaces  before  they  arrive  at  the  cryopump. 
One  seeks  a  configuration  which  can  be  cheaply  fabricated  but  has  good 
efficiency  in  effecting  energy  loss  from  the  molecules  while  passing  on 
enough  molecules  for  the  cryopumps  to  function  at  reasonable  speeds 
without  overheating. 

Of  course,  the  real  baffles  are  three-dimensional,  but  they  are 
generally  of  plane  strip  elements  so  that  except  for  their  finite  length 
and  their  ends,  they  may  be  treated  in  section  (i.  e. ,  as  two-dimensional 
systems)  provided  their  length  is  large  relative  to  their  other  dimensions. 
To  precisely  account  for  end  surfaces  requires  a  three-dimensional 
approach.  Since  results  are  available  in  the  literature  for  Monte  Carlo 
runs  for  two-dimensional  baffle  geometry,  the  M.  K.  system  can  be 
compared  to  Monte  Carlo  results  for  specific  cases. 

We  shall  go  into  some  detail  for  the  configuration  of  slanted  parallel 
plates  known  as  "louvres"  and  their  back-to-back  arrangement  called 
"chevrons"  and  examine  their  behavior  by  the  molecular  kinetics  method. 

The  louvres  considered  are  optically  tight  for  normal  incidence; 
the  chevrons  are  opaque  for  any  incidence  angle  (Fig.  1). 

For  the  detailed  study  of  the  louvre,  all  faces  (entrance  to  left,  exit 
to  right)  were  divided  into  eight  equal  zones,  and  the  resulting  32  by  32 
probability  table  was  computed  (Table  1).  We  obtained  the  asymptotic 
hit  vector  for  input  on  each  separate  zone  of  the  entrance  face  when  the 
entrance  ard  exit  faces  have  zero  reflection  coefficient  and  all  other 
faces  have  unity  (perfectly  absorbing  ends,  perfectly  diffusing- reflecting 
walls). 

It  is  plain  that  the  reflectivity  of  the  configuration  is  the  ratio  of  the 
total  hits  absorbed  by  the  entrance  face  to  the  total  input.  Also  the  net 
transmissivity  is  the  total  of  hits  absorbed  on  the  exit  face.  Recall  that 
total  hits  are  available  for  each  zone.  We  may  then  plot  the  reflectivity 
and  transmissivity  of  the  configuration  as  a  function  of  position  on  the  in¬ 
put  face.  Our  specific  case  has  6  =  45  deg. 

Figure  2  shows  the  distribution  of  reflected  particles  for  each  single 
input  zone;  Fig.  3  shows  the  transmission  distribution.  These  are  quite 
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sensitive  to  location  of  input;  consequently,  we  deduce  that  the  Monte 
Carlo  solution  to  this  problem  would  be  very  sensitive  to  any  bias  of  in¬ 
put  points.  For  completeness,  we  show  the  total  hit  distribution  on  the 
louvre  walls  (Fig.  4a  and  b)  which  correspond  to  pressure  distributions. 

These  results  are  for  perfectly  random  flow  incident  from  the  left. 

Since  by  means  of  the  original  distribution  vector  we  can  place  the  initial 
distribution  wherever  we  please,  and  we  examine  the  results  for  directed 
streams.  Part  or  all  of  such  streams  are  started  on  the  louvre  walls. 

The  dependence  of  reflection  and  transmission  on  stream  direction  is 
shown  in  Fig.  5,  distributed  over  the  output  (entrance)  faces  for  reflection. 

Figure  6a  shows  the  reflection-transmission  properties  for  the  45-deg 
louvre  as  a  whole  and  as  a  function  of  input  direction.  In  Fig.  6b,  this 
has  been  properly  weighted  to  account  for  the  finite  entrance  (cosine 
biased). 

Finally,  Fig.  7  shows  the  net  reflection-transmission  properties  for 
the  louvre  as  a  function  of  input  position  in  random  flow  with  the  final, 
lumped  values  marked  under  uniformly  random  input.  This  shows  that  for 
this  configuration  a  ''flat”  distribution  of  locations  on  the  entrance  face  in 
a  Monte  Carlo  run  could  be  replaced  by  a  single  "mean"  entry  point,  almost 
at  the  middle  of  the  entrance  face.  For  the  Monte  Carlo  method,  this 
would  remove  one  random  variable  from  the  process,  thereby  speeding  the 
convergence  of  the  process. 

Figure  8  is  a  summary  of  results  for  several  louvre  angles  and  a 
parallel  plate. 

3 .2.1.2  Two-Dinwnsional  Transfer  Probabilities 

We  define  a  source  zone  as  a  line  segment  extending  from  yj  =  0  to 
yj  =  A  at  xj  =  0.  It  has  an  aperture  height  A.  We  will  find  the  probability 
for  perfectly  random -oriented  and  positioned  paths  from  the  left  of  the 
aperture  (  x  <  0)  to  pass  through  the  aperture  and  above  a  point  (x2,  y2)  to 
its  right  (Fig.  9). 

We  evaluate  the  integral  for  the  fraction  passing  above  the  point  (x2,  y2) 

ta.<!)  -  -Jr  y  r~f  ds,~u,  n> 

1  A  Jt,  =  0  [_  J(2)  r»»  J 

Since  we  require  only  the  region  above  the  point  (x2y2)  we  may,  for 
the  inner  integral,  substitute  integration  along  a  circular  arc  from^mj[n  to 
W  2,  on  which  cos  <f>  2  =  1  and  dS2  /  r  j  2  =  d  <£  1 . 
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TA.m  '  tx  JT[  t.,„  j  d>'* 


*  X 


(8) 


(1  -  sin  <£min)  dy» 


By  construction  sin  ^mi,,  = 
hence 

Ta 


y»  -  yi 


(9) 


V*,1  +  ( y2  -  y,  )* 

.(2)  =  ~2x[y‘  _  *•  +  (y’  ~  y,)’  J  _0 

=  -yj-  [ A  “  V  X,*  +  (yj  -  A)*  +  \/TJ  +  y,1  J 
In  particular  if  A  is  taken  as  unity  (if  x,  y  are  scaled  to  the  unit  A) 


Tl  .(2)  =  -y-  -  V  *1*  +  (yi  -  l)1  +  v  *i’  +  y»*  J 


(10) 


By  noting  that  the  radicals  are  simply  distances  between  (*2^2^  anc* 
respectively,  (0,  A)  and  (0,  0),  we  can  use  distances  between  points 


Ta. (5) 

(See  Fig.  9c. ) 


J_  T  f0A  ~  rA2  +  r02 

2  L  '0A  J 


(ID 


Since  is  the  aperture  A,  constant,  we  can  establish  that  lines  of 
constant  (£)  occur  for  points  (x2  y2)  satisfying 

r02  “  rA2  “  constant  (12) 

and  such  lines  of  constant  /£)  are  a  family  of  confocal  hyperbolae, 
with  foci  at  (0,  0)  and  (0,  A).  ’  A  graph  of  constant  T  lines  has  some 
utility,  as  will  be  shown. 


Suppose  now  that  we  introduce  another  point  (x3yg)  and  calculate 


Ta.  (3) 


1  T'OA  -  rA2  +  r02 

2  L  ^  J 


(13) 


Provided  only  that  the  entire  new  line  segment  ^3  is  visible  from  all 
of  tqa*  it  is  plain  that  the  probability  of  random  rays  through  tqa  striking 
the  line  r23  is  just 

PoA.23  -  |  Ta,(3)  -  Ta  ,  (?)  I 

(14) 

1  I  ifA2  ~  f02  ^  ~  ^ f  A3  ~  r03  ^  I 
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where  the  absolute  value  is  indicated  to  ensure  a  positive  (meaningful) 
result,  independent  of  which  direction  on  r23  is  taken  to  be  positive. 

This  formulation,  incidentally,  is  now  independent  of  position  and 
orientation  of  the  figure  since  only  displacements  between  points  are  in¬ 
volved. 

The  confocal  hyperbola  graph  can  now  be  used  to  evaluate  transfer 
probabilities  between  two  line  elements  by  constructing  the  elements  to 
scale.  A  single  graph  for  A  =  1  is  sufficient  provided  all  distances  are 
converted  to  units  of  A  (Fig.  10). 

It  is  now  possible  to  calculate  transfer  probabilities  for  any  configu¬ 
ration  of  line  elements  in  the  plane  so  long  as  they  are  all  completely 
concave  to  one  another.  That  is,  each  one  must  be  entirely  visible  from 
all  points  of  any  other  element.  This  condition  is  obviously  satisfied  for 
any  parallel  line  segments,  hence  for  the  louvre  configuration. 

However,  it  is  necessarily  assumed  in  this  development  that  a  source 
segment  is  perfectly  randomly  filled  with  transfer  rays,  hence,  uniformly 
filled.  It  is  also  plain  that,  after  transfer  of  molecules  from  source  to 
target  arrays,  the  distributions  on  the  targets  will  be  non-uniform.  This 
dichotomy  can  be  resolved,  in  part,  only  by  taking  sufficiently  small  line 
segments,  even  for  linear  arrays.  By  this  means,  the  discrepancy  in 
assumed  distributions  and  accurate  (continuously  varying)  distributions 
may  be  made  as  small  as  we  please. 

Referring  again  to  the  loci  of  constant  T^  t™,  some  of  their  prop¬ 
erties  are  of  interest.  The  slope  of  the  hyperbola  at  point  (x2y£)  can  be 
easily  found  to  be  that  of  the  bisector  of  the  angle  included  between  the 
lines  tq2  and  rA2-  As  *be  Pobit  (x2y2)  is  removed  to  remote  distances, 
such  that  tqa  becomes  of  differential  size  relative  to  the  hyperbola 
of  constant  T^#  (2)  approaches  asymptotically  the  line  from  the  midpoint 
of  aperture  A,  inclined  at  angle  <f>m,  appearing  in 

-  -T  [  »  -  -i"  +m]  (15) 

This  relation  is  useful  in  estimates  of  probabilities  at  large  distances. 

For  transfer  to  a  line  segment  r23  at  great  distances 

pA,(f„)w  -  ■—  I  •'»  4>t  -  •»  4>t  I  (16) 


a  form  which  is  well  known  as  the  plane  "angle  factor"  from  a  differential 
line  segment  to  any  arc  between  angles  <t>2  and  ^3  or  for  a  differential  area 
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to  a  general  cylinder  whose  generating  lines  are  parallel  to  the  differen¬ 
tial  area  bounded  by  angles  •t>2>  ^3  relative  to  the  normal  to  the  differen¬ 
tial  area. 

3.2.1 .3  Transfer  Probabilities  with  Masking 

The  louvre  probabilities  are  quite  simply  evaluated  since  all  line 
elements  are  strictly  concave  relative  to  one  another.  When  there  exists 
convexity  in  the  figure,  it  is  necessary  to  modify  the  procedure  to  some 
extent.  We  examine  the  meaning  of  the  terms  convexity  and  concavity. 

Suppose  line  elements  (A,  B)  and  (C,  D)  are  so  oriented  that  (A,  B) 
extended  intersects  between  (C,  D).  (See  Fig.  11. )  The  plainly  the 
two  "sides"  of  a  (A,  B)  "see"  separate  parts  of  (C,  D).  We  must  then 
clearly  indicate  which  faces  are  the  front  faces  and  assign  zero  probabil¬ 
ities  for  relative  orientations  which  involve  back  faces.  The  point  of 
intersection  is  labeled  E,  and  E  is  an  end  point  for  a  smaller  segment 
which  may  be  used  to  compute  the  probability  (see  Fig.  14). 

This  is  the  simplest  form  of  masking  where  one  line  element  effec¬ 
tively  casts  a  shadow  on  another,  reducing  its  transmission  to  the  other. 

If  the  intersection  E  lies,  on  the  source  element,  a  fraction  f  of  the  source 
sees  the  front  face  of  the  target.  By  the  assumed  uniform  distribution 
within  an  element,  f  enters  as  a  multiplier  of  the  probability  computed  for 
the  fractional  part  of  the  source  that  does  not  see  the  target.  However,  if 
the  point  E  lies  in  the  target,  the  probability  is  computed  from  source  to 
that  part  of  the  target  in  view  from  the  front  of  the  source.  These  two 
quantities  differ  in  that  for  E  in  the  source,  only  part  of  the  source 
"radiation"  can  be  transmitted,  whereas  with  E  in  the  target,  all  of  the 
source  radiation  is  transmissible  to  part  of  the  target. 

A  more  complex  form  of  masking  occurs  when  a  convex  part  (Fig.  12) 
of  the  boundary  is  interposed  between  source  and  target.  Let  the  end  of 
the  mask  point  be  labeled  M.  Then  lines  from  A,  B,  C,  or  D  through  M 
will  cross  within  the  elements  (A,  B)  or  (C,  D)  if  masking  occurs. 

Consider  the  configuration  shown  where  (A',  B')  is  the  segment  of 
(A,  B)  where  masking  occurs.  We  see  that  (A,  A')  does  not  see  (C,  D)  at 
all,  but  that  (B\  B)  sees  all  of  (C,  D).  Thus,  (B',  B)  and  (C,  D)  may  be 
simply  inter- related  as  a  simple  concave  pair. 

It  is  plain  that  any  rays  drawn  from  (A1,  B')  to  (C,  D)  must  pass 
through  the  line  (M,  D).  Thus  we  may  use  (A',  B')  and  (M,  D)  as  a  simple 
concave  pair  to  obtain  the  probability  from  (A',  B')  to  (C,  D).  Likewise, 
all  rays  from  (C,  D)  to  (A,  B')  must  pass  through  (M,  B').  So  the  probability 
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for  (C,  D)  to  (A,  B')  is  computed  from  the  simple  concave  pair,  (C,  D) 
to  (M,  B')-  To  this  is  added  the  contribution  of  (C,  D)  to  (B,  B').  This 
situation  is  summarized  in  Fig.  12. 

More  complicated  masking  can  occur  in  some  configurations  such 
as  two  or  more  points  like  M.  They  may  all  be  handled  by  arguments 
similar  to  these,  breaking  the  segments  into  more  simply  related  parts. 

We  are  how  in  position  to  compute  all  transfer  probabilities  for 
the  chevron  configuration,  as  illustrated  in  Fig.  13.  We  use  Roman 
numerals  to  distinguish  the  six  major  faces  and  use  a  symmetrical  num¬ 
bering  scheme  for  the  six  N  faces,  N  being  the  number  of  segments  in 
each  face.  (It  is,  of  course,  not  necessary  to  use  either  equal  size 
segments  or  equal  numbers  of  segments  per  face. ) 

We  see  that  probabilities  between  faces  I,  II,  and  TII  are  already 
obtained  in  the  louvre  case.  (By  symmetry,  whatever  we  say  of  faces 
I,  II,  and  III  holds  for  IV,  V,  and  VI. )  Also  the  upper  half  of  face  IV 
is  simply  concave  toward  I,  II,  and  in,  while  V  and  VI  are  not  visible 
from  I  and  II.  Part  of  III  sees  IV,  V,  and  VI,  but  the  relations  between 
I,  IV,  and  VI  (for  lower  in)  involve  masking  by  the  junction  of  II  and  V. 

We  note  also  that  the  first  type  of  masking  discussed  can  occur  for  II, 

IV  transfers,  but  this  is  removable  by  making  the  point  of  intersection 
of  II  (extended)  with  IV,  a  point  of  division  segments  of  IV.  For 
6  =  45  deg,  we  provide  this  by  requiring  an  even  number  of  segments 
per  face. 

3.2.1  A  Criticism  of  Louvro  and  Chovron  Results 

By  comparison  with  results  of  other  workers  using  Monte  Carlo  or 
laboratory  methods,  the  present  transmission  factors  are  usually  a  few 
percent  high.  This  is  a  direct  consequence  of  the  assumed  uniformity 
of  distribution  of  flux  on  the  target  zone,  an  assumption  which  is  in  reality 
violated  at  any  concave  corner.  In  effect,  the  original  flux  is  inserted 
somewhat  deeper  into  the  figure  than  intended. 

The  following  illustrates  the  point.  Suppose  we  had  two  adjacent 
segments  at  right  angles.  We  find  P  =  0.  283  as  the  transfer  probability. 
Now  reduce  segment  size  by  one  half,  and  the  probability  is  unchanged. 
Thus,  at  the  corner,  if  the  segment  is  1/n  of  the  original  segment,  then 
0.  283/n  of  the  original  flux  is  transferred  at  the  corner,  and  less,  else¬ 
where.  The  corner  is  a  "hot  spot"  with  highest  concentration.  Now  the 
assumption  of  uniformity  in  a  target  zone  plainly  causes  an  effective  shift 
of  the  flux  away  from  the  comer. 
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This  error,  commonly  found  in  radiant  heat-transfer  numerical 
calculations,  can  be  reduced  only  by  finer  subdivision,  particularly  at 
corners  but  generally  in  regions  where  the  transfer  probability  changes 
most  rapidly  with  location.  By  this  means,  not  only  are  the  distribu¬ 
tions  shifted  a  lesser  amount,  but  the  number  of  zones  causing  signifi¬ 
cant  error  is  relatively  reduced. 

3.2.2  Thraa-Dimantioiial  Application* 

322.1  Gaomatry 

The  problem  of  concentric  spheres  constituted  our  introduction  to 
the  M.  K.  method,  a  problem  specially  formulated  for  its  application 
from  which  we  have  generalized  the  method  to  other  configurations. 

The  investigation  involved  developing  a  concept  of  "pressure"  in 
molecular  flow,  in  particular  for  a  spherical  space  environment  simula¬ 
tor.  An  inner  concentric  sphere  was  an  idealized  space  vehicle.  The 
gas  flow  space  was  thus  completely  bounded  between  the  outer  "chamber" 
and  the  inner  "vehicle". 

The  problem  was  formulated  by  C.  Tsonis  of  General  Electric  Com¬ 
pany  (under  contract  AF  40(600)-954)  with  computing  support  provided  by 
ARO,  Inc. ,  Arnold  Engineering  Development  Center.  *  Tsonis  used  the 
distinction  between  the  concave  chamber  and  the  convex  vehicle  to  sepa¬ 
rate  the  transfer  probabilities  into  three  parts  (effectively  partitioning 
the  probability  matrix)  and  invoked  the  symmetry  of  the  system  to  reduce 
the  calculating  effort.  Most  significant  of  all,  Tsonis  integrated  the 
integral  expressions  for  transfer  probabilities  to  provide  closed  forms 
for  direct  numerical  calculations. 

As  originally  proposed,  our  problem  was  simply  to  program  the  ex¬ 
pressions  provided.  It  was  intended  that  the  spheres  be  divided  into  2N 
coaxial  spherical  zones,  where  N  was  a  "large"  number,  nominally  20  to 
40.  Since  three  kinds  of  transfer  occur,  this  means  we  had  to  provide 
storage  for  three  80  by  80  arrays,  in  effect.  Equatorial  symmetry  of 
geometry  reduced  this  to  three  40  by  80  arrays,  and  a  simple  relation  be¬ 
tween  two  of  these  resulted  in  only  two  40  by  80  arrays  being  required. 
This  made  the  program  tractable  for  use  with  our  computer  (IBM  7070, 
10K  storage).  The  program  was  blocked  to  compute  probabilities  first, 
then  to  compute  the  flow  characteristics  with  another  program. 


*C.  A.  Tsonis,  "Molecular  Flux  Distribution  in  an  Aerospace  Chamber 
Analysis  of  Gas  Kinetics  -  Summary  Report,  "  AEDC-TDR-63-88, 

April  1963. 
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The  "bounce"  technique  was  first  utilized  with  control  option  for 
the  number  of  bounces.  Later,  the  same  problem  was  proposed  to  be 
done  in  terms  of  matrix  algebra  to  arrive  directly  at  asymptotic  values. 

For  the  matrix  method,  the  requirement  of  computing  speed  and  the 
necessity  for  more  arrays  in  storage  led  to  reducing  the  program  capa¬ 
bility  so  that  the  maximum  array  was  40  by  40,  which  proved  adequate. 

In  every  case  so  far  tried,  the  matrix  (40  by  40  maximum)  to  be  inverted 
was  so  well  conditioned  that  single  precision  matrix  inversion  (by 
diagonalizing)  was  adequate. 

The  equations  programmed  are  included  in  this  paper  (Fig.  14)  to¬ 
gether  with  a  special  demonstration  of  pulse  propagation  through  a 
small  space  between  a  chamber  and  a  vehicle  that  nearly  fills  the  cham¬ 
ber  (Fig.  15). 

For  plane  surface  boundaries  and  for  some  simple  shapes  in  various 
relative  positions,  the  transfer  probabilities  may  be  calculated  from  the 
formulae  in  the  literature.  More  generally,  one  can  expect  the  neces¬ 
sity  of  developing  numerical  values  from  the  fundamental  integral  definition. 

D  _  1  f  f  cos  4>x  cos  <f>2  dAA,  dAA, 

aa.aa,  aa,  J\Al  J\Kt  n 

It  is  quite  possible  to  numerically  integrate  this  quadruple  integral; 
the  integration  limits  often  cause  considerable  difficulty. 

It  should  be  noted  that  integration  can  also  be  done  by  the  Monte 
Carlo  method  by  estimating  the  fraction  of  starts  from  the  source  that  will 
strike  the  target  segment  somewhere.  This  is  quite  feasible  since  only 
single  step  processes  are  involved,  not  long  process  chains. 

We  will  devote  attention  to  the  concentric  sphere  case.  For  axisym- 
metric  geometry  with  constant  latitude  zone  boundaries,  the  probability 
integrals  were  evaluated  by  Tsonis  and  are  given  in  Fig.  14. 

3 .2.2.2  Transfer  Equations 

The  relative  size  of  two  concentric  spheres  is  completely  defined  by 
the  relation  rinner  /  router  =  sin  y.  We  chose  a  polar  coordinate  system 
with  a  common  polar  axis.  Let  each  sphere  be  divided  into  2N  zones 
whose  boundaries  are  plane  circles  parallel  to  the  equatorial  plane,  de¬ 
fined  by  constant  angle  0  from  the  equator  and  having  uniform  angular 
width  A  0  =  x/2N. 

A  zone  is  associated  with  its  lower  boundary,  and  its  zone  index  num¬ 
bers  (I  for  source  zone,  J  for  target  zone)  start  at  the  south  pole,  0  =  -  jr/2. 
Theta  (0)  is  the  latitude  angle. 
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We  assume  axial  symmetry  so  that  our  expressions  hold  uniformly 
in  longitude.  Hence  the  probability  of  a  molecular  transfer  from 
zone  I  to  zone  J  depends  solely  on  0j,  0j,  Ad,  and  y. 


In  obtaining  the  transfer  probabilities,  transfers  from  zone  0i  to 
the  entire  region  above  0j  are  calculated: 

P  («i,  e  ;>  0}) 


and 

P  (*i.  0  £  0j  +  i ) 


(Fig.  16). 
then 

P  ( 0, .  0j)  =  P  (0i.  0  £  0j)  -  P  (0j,  0  £  0j  +  l ) 


for  all  0j  but  for  only  0  £  9j  <  irf  2.  By  equatorial  symmetry  the  remain¬ 
ing  values  are  assigned  since  the  symmetry  means  that 

P (01  0j)  -  P (02N  -1+1,  02N  -  J  +  l) 


when  I  and  J  range  from  1  to  2N. 

There  are  three  types  of  transfers  (Fig.  17):  outer  sphere  to  inner 
sphere  (1),  inner  sphere  to  outer  sphere  (2),  and  outer  sphere  to  outer 
sphere  (3),  The  type  inner  to  inner  does  not  occur  since  the  inner  sphere 
is  totally  convex  to  itself.  In  all  types  of  transfers  we  limit  the  calculation 
to  transfers  from  any  source  location  to  "caps"  above  the  equator  by 
reason  of  equatorial  symmetry. 

In  the  outer  sphere  to  inner  sphere  transfers  (1),  no  part  of  the  upper 
hemisphere  is  visible  from  latitudes  below  0j  =  x/2  -  y.  The  part  of  the 
inner  sphere  visible  from  a  point  at  0j  >  jt/2  -  y  is  the  region  in  common 
with  the  cap  (0j)  and  the  cone  of  half  vertex  ang*e  y.  This  is  described  by 
Eq.  (2)  of  Fig.  14  until  0j  >  y.  When  this  point  occurs,  some  caps  (0j)  lie 
wholly  within  the  y  cone  of  the  point  at  0j,  and  these  wholly  contained  caps 
have  probabilities  calculated  from  Eq.  (1)  of  Fig.  14  (Fig.  18a). 

In  the  inner  sphere  to  outer  sphere  transfers  (2),  the  bounding  curves 
are  the  cap  (0j)  lower  boundaries  and  the  horizon  circle  on  the  outer  sphere, 
cut  by  the  plane  that  is  tangent  to  the  point  (fij)  on  the  inner  sphere.  The 
same  limits  on  0j  apply  as  for  type  (2)  transfers.  Near  the  south  pole, 
there  is  a  range  of  0j  for  which  the  horizon  circle  does  not  project  above  the 
equator.  Then  there  is  a  range  such  that  the  horizon  circle  cuts  caps  in 
the  northern  hemisphere  so  that  the  cap  boundary  and  the  horizon  circle  to¬ 
gether  define  the  region  to  which  transfer  may  occur,  and  Eq.  (4)  of 
Fig.  14  is  used.  Finally,  some  caps  occur  wholly  contained  within  the 


15 


AEDC.TDR-63-120 


horizon  circle,  and  these  have  probabilities  calculated  from  Eq.  (3)  of 
(Figs.  18b  and  14). 

In  the  outer  sphere  to  outer  sphere  transfers  (3),  the  main  boundary 
feature  is  the  shadow  of  the  inner  sphere  cast  on  the  outer  sphere 
(Figs.  18c  and  14). 

If  the  shadow  is  wholly  below  the  spherical  cap  (0j),  the  probability 
is  given  by  Eq.  (5).  If  this  shadow  cuts  the  spherical  cap  (0j)  then 
Eq.  (7)  is  used.  But  if  the  shadow  is  wholly  contained  within  the  cap 
(6j) ,  Eq.  (6)  is  used. 

In  the  cases  computed,  y  was  taken  to  be  some  multiple,  N2,  of  A 0, 
so  that  all  tests  concerning  which  equation  to  use  are  based  on  relations 
between  zone  index  numbers  I  and  J  and  the  numbers  N2  and  N.  It  is 
well  to  note  that  Eqs.  (1)  and  (2)  differ  from  Eqs.  (3)  and  (4)  only  by  a 
factor  sin^  y ;  the  computer  programs  take  advantage  of  this  as  well  as 
symmetry  to  minimize  storage  problems. 

3.2.2.3  Example  of  Results  Obtained 


We  show  (Fig.  15)  a  plot  of  results  using  the  bounce  system  to  de¬ 
scribe  the  propagation  of  a  pulse  of  molecules  originating  on  the  zones 
on  either  side  of  the  equator.  The  same  input  is  started  on  both  spheres, 
and  the  result  is  a  wave  propagating  toward  the  poles.  A  reflection  coef¬ 
ficient  of  unity  preserves  all  the  molecular  flux.  It  is  apparent  that  the 
pulse  is  diffusing  as  well  as  moving  away;  after  several  hundred  bounces 
the  gas  would  be  practically  uniformly  distributed. 

For  further  applications  see  Tsonis  (AEDC-TDR- 63-88)- 
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APPENDIX 


THE  "COSINE  BIAS"  LAW 

Let  us  imagine  a  space  in  which  particles  have  completely  random 
positions  and  directions  of  motion.  If  we  now  introduce  a  plane  imagi¬ 
nary  surface  into  this  space  and  examine  the  passage  of  particles 
through  the  surface,  we  see  the  controlling  factor  is  the  projected  area 
of  the  surface.  The  maximum  number  o*  paths  through  the  surface  is 
associated  with  directions  normal  to  the  surface,  and  minimum,  for 
directions  parallel  to  it.  The  projection  of  an  area  in  a  given  direction 
invokes  the  cosine  of  the  angle  between  that  direction  and  the  normal  to 
the  area.  This  is  the  purely  geometrical  origin  of  the  cosine  bias  law. 

To  see  that  this  law  holds  in  two-dimensional  geometry,  it  suffices 
to  observe  that  completely  random  directions  in  three  dimensions  pro¬ 
ject  onto  any  plane  as  equally  random  directions  and  that  the  projection 
of  a  line  segment  in  a  given  direction  still  involves  a  cosine  factor. 

In  either  case,  in  writing  a  cosine  bias  distribution  for  directional 
probabilities,  it  is  mandatory  to  normalize  the  distribution  so  that  the 
integral  of  the  distribution  over  all  realizable  directions  is  unity. 

We  note  that  the  above  discussion  is  independent  of  rates  of  motion, 
i.  e. ,  of  velocity,  and  the  cosine  bias  law  does  not  incorporate  any  in¬ 
crease  of  "forward"  scattering  at  very  oblique  incidence. 


ALTERNATE  PLANE  PROBABILITY  CALCULATION 

We  will  develop  here  an  alternative  approach  to  calculating  transfer 
probabilities  in  the  plane.  We  will  find  the  average  fraction  of  parallel 
streams  from  a  source  aperture  to  the  region  above  a  given  point,  averag¬ 
ing  over  all  orientations  of  the  stream. 

As  before,  the  source  is  a  line  segment  on  the  y-axis,  0  <  y  <  A.  A 
point  (x,  y)  is  to  the  right.  The  source  is  illuminated  from  the  left  with 
a  strictly  parallel  beam  inclined  at  angle  <t>  to  the  x-axis.  The  beam  at 
least  fills  the  source  aperture  (Fig.  19). 

The  point  (x,  y)  may  be  in  or  out  of  the  beam  issuing  from  the  source. 
We  calculate  the  fraction  f {<f>,  x,  y)  of  the  beam  passing  above  (x,  y).  The 
width  of  beam  through  (0,  A)  is  A  cos  <f>. 
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Plainly,  for  general  finite  (x,  y)  anywhere  to  the  right  of  (0,  A), 
there  is  some  range  of  <£  for  which  no  part  of  the  beam  passes  above 
(x,  y);  this  is  the  range  -  (nl 2)  <  <t>  <  (x,  y.  A).  Also  there  is 

some  range  for  which  the  entire  beam  passes  above  (x,  y);  this  is  the 
range  <jf>max  (x,  y,  0)  <  <  (it 1 2).  And  finally,  over  the  range 

«£min<  ^  <  ^max«  some  fraction  of  the  beam  passes  above  (x,  y). 

We  define  the  two  distances 

r«(x,y,A)  -  Vx1  +  (y  -  A)’ 


r0(x,y,o)  -  y/x1  +  y* 

and  from  the  figure,  find  the  beam  width 

A  cos  <f>  =  rQ  sin  (<£max  -  <f>)  +  rA  sin  (<j>  -  <^min) 

and  the  fraction  of  beam  passing  above  (x,  y)  is  thus 
f  -  0 


■N  f-  -y  <  0  <  <t> min 


!  -  (r  .  sin  ( <j>  ~  <f> min  )  )  /  A  coe  <f> 


f  -  i 


<£min  <  <f>  <  <t>  max 
<t> max  <  $  <  -§" 


We  wish  to  find  the  average  fraction  f  of  a  beam  of  width  A  cos  4>  as  4> 
ranges  from  -  (jr/2)  to  (jt/2),  given  by  the  expression 


rr 

T- 


f  = 


J  n  A  f  (<i>)  cos  <f>  <•  4> 

t'—T _ 

n 

T 


|  A  cos  <f>  d  <f> 

*—  ( ...  <J>  Z^Aa  co.  4  J  *  .  f 


1 

2  A 


~r 


-*r- 


r.  cos  ( <Amtx  -  <£min)  +  r.  +  A  -  A  sin  <t>  m%x 


■] 


From  the  definitions  of  r^  and  rQ,  we  find  x/ ra  =  cos  #min» 
x/rD  =  cos  *max.  y/rQ  =  sin*max,  (y- A)/ta  =  sin  *min,  so 
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f  -  |^A  ( 1  -  cos  <t>max  cos  <£min  -  sin  <£max  sin 

+  A  (1  -  sin  <£raax)J 

■  -Sir1)  **(■-*)] 

♦*] 

■  Ta[U  '  *  '*] 

This  is  the  form  derived  previously  by  a  quite  different  approach. 


THE  CHARACTER  OF  THE  PROBABILITY  MATRICES 

For  the  linear  two-dimensional  arrays,  louvre  and  chevron  (Fig.  1), 
it  is  plain  that  the  equations  for  Py  give  zero  for  collinear  segment 
probabilities.  The  matrix  will  thus  always  be  zero  on  the  principal  diag¬ 
onal  and  for  groups  (minors)  in  any  one  face.  Other  zeroes  may  occur  for 
complete  masking. 

There  are  also  in  these  configurations  symmetrical  properties  which 
appear  in  the  matrix  if  only  a  symmetrical  numbering  scheme  is  used. 

We  saw  previously  that  the  matrix  to  be  inverted  (for  asymptotic 
values)  is  of  form 


I  I  -  A  P ' I  or  I  I  -  P'A  I 

which,  due  to  the  collinearity  mentioned  above,  is  always  unity  for  every 
main  diagonal  element.  Since  we  have  °-Ay  v  and  0  <  Pjj  <  I  in  all 
open  polygonal  configurations,  the  matrix  to  be  inverted  always  has  a 
strong  "backbone".  This  explains  why  matrix  inversion  works  very  well 
even  for  48  by  48  matrices  and  single  precision  in  these  problems. 

In  more  general,  say  three-dimensional  cases,  the  diagonal  (Py) 
elements  will  be  non- zero  if  a  surface  element  is  actually  concave  (as 
within  a  segment  of  a  sphere).  By  further  subdivision,  they  may  be  made 
quite  small,  while  the  (I  -  AP')  or  (I  -  P'A)  matrix  maintains  a  strong 
trace.  Thus  the  only  trouble  expected  from  the  matrix  inversion  would 
occur  for  extremely  numerous  divisions;  this  is  round-off  error,  and 
even  this  limitation  can  be  extended  by  using  double  precision  matrix  in¬ 
version. 
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1  2  3  4  5  6  7  8 

Output  Zone 

Ft*.  2  Raflactian  Diatrlbutiau  by  Zona*  far  Sia*la  Zana  Input  (Lauvra) 
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Output  Zona 

Fi|.  3  TmtalulM  DlitriWliM  by  Zona*  far  Slnpla  Zona  Input  (Lauvra) 
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b.  Top  Wall 
Fif.  4  Cooclodod 
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8  o.  x 


-t«n~1(l/7)  -  |  -  <t> 

— tan"1 (2/6) 

- - 7 

r—  tan-1 (3/S) 


-tan"1  (4/4) 


.tan"1 (5/3) 


tan"1  (7/1) 
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Fig.  7  Total  Output  for  Singlo  Input  Zono  (Louvre),  0  -  45  dog 
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Fig.  10  Curves  of  Constant  Transmission 
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A 


P(AB  •*  CD)  -  P(AB  -»  ED) 
FI®.  11  SimpU  Masking 
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P(AB  «*■  CD) 


.A7!!  x  P(A'B' 

XB 


-*  MD)  + 


x 


pCePb  -*  CD) 


B  D 


P(CB  •*  AB) 


_  p(CD  -*  B  B)  +  P(CD  MB  ) 


Fig.  12  Ma»king  by  Cmrm*  Con»t 
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8  -  Ancl*  defining  loner  boundary  of  polar  none 
y  -  Arc  ala  *^/*c 

I  -  Index  of  nource  none 
J  -  Index  of  target  zone 
Subecrlpt  c  -  "chaaber  wall",  outer  epbere 
v  -  “vehicle",  Inner  apbere 


1. 


a. 


*cv<W  - 


nln  Sj  a In  y  „T 
- , - T 


V 


nln  6j  nln  7 
- 27 - 


*  -  CO 


-ir,ia  y  ~ b1b  fli b1b  Vi 

I  coe  e,  coe  Uj  I 

,r«ln  8.  -  nln  y  ala  e.-i 
®o  "  BlB  t  coe  y  coe  &l  J 

•In  y  f 1  ♦  nln*  y)  nln  9j  -  2  nln  y  nln  flj 
»|  [a  nln  y  con(0j  ♦  8j)  ♦  1  ♦  nln*  ij  [-1  nln  y  con  (8j  -  #j)  ♦  1  ♦  nln^  y] 

peon  <8X  -  8j)  -  nln  Y[U* 

X1  “Icon  l8j  -  +  nln  yj 

x_x  3  nln  y  com  <8j  ♦  8j)  ♦  1  ♦  nln*  y 

***  ^^lZ^~com~(e^~e^)~*T7~mii?~y 


•la  8, 


3-  'nc<‘i*V  “  rar?  “  — !S- 


'  2  •ina  y 


<•  P,c(W  "  It  hiny  *0  *  *  Po)  '  tBB_1  [Xi(£t)  ] 

<*0.  P0-  °>  xl*  T?T  *boT*) 


5.  Pee(«,.«j) 


*•  Pcc<9I'V  - 


•In  9 


♦  a 


•la  9, 


•  In  Y 


2.  P. 


.<  w  -  -  ^  t  . 

jf  coe  3y  ♦  ala  9{  nln  82 
*0  ‘  COB’  [*  - con  Bj  coe  S} - j 

-lT  BiB  ®J  1 

p0  -  eln  Icot  ay  tan  8j  ♦  ,U  Ty  coe  9,J 

[coe  (9j  -  8j)  ♦  coe  ay] 1/1 
X1  ”  |co«  l8j  ♦  Vj>  -  con  ayj 


(&t) 


1/a  1  -  con  (8j  ♦  9j> 

’  am  «j  ♦  eia  VJ 


Fig.  14  Probability  Equation*  of  Concentric  Sphere* 
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Fig.  15  Puls*  Propagation  botwoon  Coneontrie  Sphoras 
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To  find  P(0^,0j),  probability  for  transfer  from  zone  i  to  zone 
j,  find  first  for  "caps" 

P (0i,  e  >  ©j) 

Then  P(6i,  6  >  0  ) 

and  P^,  ej)  “  P(ei.  0  >  0j)  -  P(0if  0  >  0j+1)  ,  for  zone  j 

For  either  P(@i,  &  >  6j)  and  PC©^  9  >  ej+1)  the  boundary  arcs 
are  and  C6j  or  C0j  +  1- 


At  left  the  "caps"  bounded  by  C0j  and 
C0j+1  are  bounded  by  these  arcs  alone, 
but  the  arcs  and  C9^  are  both 

involved  for  lower  zones. 


Fig.  16  Zone  CtOMtry,  Concentric  Spheres 
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Chamber  Wall 
to  Vehicle 


Chamber  Wall  to 
Chamber  Wall 

Shadow  of 
Vehicle 


Fig.  17  Melecwler  Transfer*  (Target  in  Northern  Hemisphere) 
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Chamber  to  Vehicle,  y  <  ^ 


N  N 


Note:  Aa  y  •*  -J  point  b  and  c 
merge;  type  III  doee  not 
occur  for  y  >  Equations 
refer  to  equations  In  Fig.  14. 


areas. 

a.  Outer  Sphere  te  laser  Sphere 
Fip.  IS  Type*  ef  Transfer* 
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Chamber  to  Chamber 


I  Zones  bounded  below  a  Include  H Zones  bounded  above  c  are 

entire  shadow,  above  b  are  free  of  shadow;  all  below 

wholly  in  shadow,  between  are  cut  shadow, 

cut  by  shadow  line. 


Equations  refer  to  equations 
in  Fig.  14. 


Fa  ( M 


HI  All  northern  hemisphere  is  S  Zones  bounded  above  e  are  free 

free  of  shadow.  of  shadow,  below  d  Include 

whole  shadow,  between  cut 
shadow . 

c.  Outer  Sphere  to  (feter  Sphere 
Fig.  18  Concluded 
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